Introducción al análsis de imagenes satelitales

Este video exploraremos cómo importar y graficar un ráster multibanda en R. Usaremos un subset de una escena de Landsat 9 producto de surferce reflectance, perteneciente al bosque Chaqueño. Primero aprenderemos a leer un raster, y visualizarrlo. Finalmente, calcularemos distintos indices de vegetación.

En este notebook vamos a utilizar el producto de reflectancia de superficie de Landsat 9. El último satélite de la NASA fue lanzado el 27 de septiembre de 2021 a bordo de Atlas V y sus imágenes están disponibles desde el 10 de febrero de 2022. Como sus sucesores, cuenta con el instrumento OLI (Operational Land Imager 2) y TIRS (Thermal Infrared Sensor 2). En este caso, OLI 2 y TIRS 2 introducen la mejora en los datos provenientes de masas de agua, humedad del suelo y vegetación.

La flota Landsat 9 continua la sistemática de bandas de las anteriores misiones. Seis bandas para el instrumento OLI 2 basadas en el visible, NIR y SWIR junto a la banda pancromática y las bandas de cirros y aerosol. Dos bandas térmicas para el TIRS 2. En este caso la resolución radiométrica pasa de 12 a 14 bits frente a los datos de Landsat 8.

A continuación, te mostramos una comparativa de las bandas de ambos satélites (Landsat 9):

Bandas

Sus resoluciones se mantienen a 30 metros en el visible, el NIR y el SWIR pudiendo trabajar la banda pancromática para un refinado que permita la mejora de la imagen hasta 15 metros de resolución (frente a los actuales 10 metros de Sentinel 2 en el visible)

El periodo de revisita de las imágenes es de aproximadamente 8 días generando tiles de aproximadamente 170 km por 183 km determinados de manera constante a través de los convencionales parámetros de Path y Row. Si necesitas conocer exactamente los momentos futuros o pasados de mapeo para Landsat 9, puedes recurrir a la herramienta online Landsat Acquisition Tool que te permitirá visualizar los mapeos temporales en un momento determinado para la flota del nuevo satélite.

Librerias a utilizar

# load the raster, sp, and rgdal packages
rm(list=ls())
library(raster)
## Loading required package: sp
library(sp)
library(rgdal)
## Please note that rgdal will be retired by the end of 2023,
## plan transition to sf/stars/terra functions using GDAL and PROJ
## at your earliest convenience.
## 
## rgdal: version: 1.5-32, (SVN revision 1176)
## Geospatial Data Abstraction Library extensions to R successfully loaded
## Loaded GDAL runtime: GDAL 3.4.2, released 2022/03/08
## Path to GDAL shared files: /Library/Frameworks/R.framework/Versions/4.2/Resources/library/rgdal/gdal
## GDAL binary built with GEOS: FALSE 
## Loaded PROJ runtime: Rel. 8.2.1, January 1st, 2022, [PJ_VERSION: 821]
## Path to PROJ shared files: /Library/Frameworks/R.framework/Versions/4.2/Resources/library/rgdal/proj
## PROJ CDN enabled: FALSE
## Linking to sp version:1.5-0
## To mute warnings of possible GDAL/OSR exportToProj4() degradation,
## use options("rgdal_show_exportToProj4_warnings"="none") before loading sp or rgdal.
library(ggplot2)
# load raster in an R object
setwd("/Users/veronica/Documents/Maestria/modelado_espacial/modelado_espacial/trabajo_final/data")

Ahora realicemos un gráfico de la imágen.

# Blue
imagen <- raster('L9_chaco.tif')
imagen_stack = stack("./L9_chaco.tif")
summary(imagen)
## Warning in .local(object, ...): summary is an estimate based on a sample of 1e+05 cells (9.86% of all cells)
##           L9_chaco
## Min.    -0.0571925
## 1st Qu.  0.0179375
## Median   0.0197800
## 3rd Qu.  0.0228050
## Max.     0.4840350
## NA's     0.0000000
# create a grayscale color palette to use for the image.
grayscale_colors <- gray.colors(100,            # number of different color levels 
                                start = 0.0,    # how black (0) to go
                                end = 1.0,      # how white (1) to go
                                gamma = 2.2,    # correction between how a digital 
                                # camera sees the world and how human eyes see it
                                alpha = NULL)   #Null=colors are not transparent

# Plot band 1

# plot all three bands separately
plot(imagen_stack, 
     col=grayscale_colors)

También podemos ver el métadato, para saber la resolución espacial, las dimenciones de lasa bandas, el nombre de las mismas y la proyección.

# view attributes: Check out dimension, CRS, resolution, values attributes, and 
# band.
imagen_stack
## class      : RasterStack 
## dimensions : 789, 1285, 1013865, 8  (nrow, ncol, ncell, nlayers)
## resolution : 0.0002694946, 0.0002694946  (x, y)
## extent     : -62.29071, -61.94441, -25.94532, -25.73269  (xmin, xmax, ymin, ymax)
## crs        : +proj=longlat +datum=WGS84 +no_defs 
## names      : SR_B1, SR_B2, SR_B3, SR_B4, SR_B5, SR_B6, SR_B7, ST_B10

Si estamos interesados en unicamente una banda, podemos importar directamente la banda de interes.

# Can specify which band we want to read in
banda_2 <- 
  raster(paste0("L9_chaco.tif"), 
           band = 2)

# plot band 2
plot(banda_2,
     col=grayscale_colors, # we already created this palette & can use it again
     axes=FALSE,
     main="RGB Imagery- SR_B2-Green")

Vamos a analizar e interpretar el contenido de información de las bandas de las imágenes ópticas multiespectrales. Familiarizarnos con las diferentes opciones de combinación de bandas y la obtención de imágenes de color.

Primero veamos como es la distribución de los valores de reflectancia de superficie de todas las bandas.

# view histogram of all 3 bands
hist(imagen_stack,
     maxpixels=ncell(imagen_stack))

No obtenemos tanta información de estos gráficos en escala de grises; a menudo se combinan en un “compuesto” para crear visualizaciones más interesantes.

Para hacer una imagen de “color verdadero (o natural)”, es decir, algo que parezca una fotografía normal (vegetación en verde, agua azul, etc.), necesitamos bandas en las regiones roja, verde y azul. Para esta imagen de Landsat, se pueden usar las bandas 4 (roja), 3 (verde) y 2 (azul). El método plotRGB se puede utilizar para combinarlos en un solo compuesto.

#landsatRGB <- stack(b4, b3, b2)
plotRGB(imagen_stack, 
        r = 'SR_B4', g = 'SR_B3', b = 'SR_B2',axes = TRUE,  main = "Landsat Composite")

Observen que la imagen de arriba se ve oscura. También puede proporcionar argumentos adicionales a plotRGB para mejorar la visualización (por ejemplo, una extensión lineal de los valores, usando strecth = “lin”).

#landsatRGB <- stack(b4, b3, b2)
plotRGB(imagen_stack, 
        r = 'SR_B4', g = 'SR_B3', b = 'SR_B2',axes = TRUE, stretch = "lin", main = "Landsat True Color Composite")

Las tonalidades más habituales en una composición en falso color son:

#landsatFCC <- stack(b5, b4, b3)
plotRGB(imagen_stack, 
        r = 'SR_B5', g = 'SR_B4', b = 'SR_B3', axes=TRUE, stretch="lin", main="Landsat False Color Composite")

## Indices de vegetación

Los índices son mediciones empíricas de las propiedades de la superficie terrestre; constituyen valores adimensionales basados en medidas de radiación y calculados a partir de la combinación de bandas espectrales.

El objetivo principal de los índices es formular una medida sintética precisa sobre las variaciones espacio-temporales que ocurren en los ecosistemas a partir de las bandas espectrales de una imagen, dado que asumen una relación (empírica o teórica) con variables bio-geofísicas de la superficie.

El caso más conocido es el Índice de Vegetación de Diferencia Normalizada (NDVI, Normalized Difference Vegetation Index) (Tucker, C. J. 1979). Las bases teóricas para su desarrollo se derivan del comportamiento de la firma espectral típica de la vegetación. La energía reflejada en el visible es muy baja debido a la actividad de los pigmentos fotosintéticos que determinan una fuerte absorción en las porciones del espectro correspondientes al azul (470 nm) y al rojo (670 nm). Por el contrario, casi toda la radiación en el infrarrojo cercano se refleja (hay muy poca absorción) dependiendo del índice de área foliar (LAI), la distribución angular de las hojas y su anatomía y morfología. El contraste entre las respuestas del rojo (R) y el infrarrojo cercano (IRc) constituye una medida sensible de la cantidad de vegetación, donde la máxima diferencia entre R y IRc corresponde al estadío de mayor densidad o vigor de la vegetación, y el mínimo contraste a áreas de muy poca vegetación o vegetación en estado senescente. Este salto se denomina “red edge” o “borde rojo”

El más sencillo de estos índices es simplemente el cociente entre las bandas roja e infrarroja:

$ VI= IRc/ R$

El NDVI, mencionado anteriormente corresponde al cociente:

$ NDVI = (nir -red) / nir +red)$

Como cociente normalizado tiene ciertas ventajas. Reduce parcialmente los efectos ambientales (atmósfera) que se observan en las bandas individuales, variaciones provocadas por la topografía, sombras, variaciones en iluminación. Realza variaciones pequeñas en las respuestas espectrales de las coberturas. Finalmente, este índice se puede determinar a partir del contaje digital crudo, radiancias, reflectancias TOA, reflectancias de superficie. Sin embargo, si bien las unidades se cancelan, y el cociente neutraliza relativamente el ruido introducido por la atmósfera, los resultados indican que es recomendable trabajar con las imágenes en reflectancia de superficie (corregidas atmosféricamente), especialmente cuando se trabaje con imágenes de más de una fecha (o con series temporales).

En algunos casos puede ser interesante reemplazar la banda del rojo visible por la del verde visible, donde entrarían en juego otros pigmentos y la cual es sensible al estrés en las plantas.

El NDVI describe adecuadamente el comportamiento de la vegetación en ambientes donde la cobertura vegetal es alta o por lo menos hay baja proporción de suelo desnudo. En casos donde la cobertura de la vegetación es muy pobre una buena parte de la reflectancia registrada por los sensores corresponde a la contribución del suelo y el NDVI no tiene un comportamiento eficiente frente a estas situaciones (Price, J. C. 1987). El Índice verde corregido por suelo (SAVI, Soil Adjusted Vegetation Index), da cuenta de este fenómeno (Huete, A. R. 1988):

$ SAVI= ((nir -red)/(nir +red+L))*(1+L)$

Donde L es la contribución del suelo y puede variar entre 0 y 1, según el porcentaje de suelo desnudo. Un valor de 0,5 indicaría una proporción de 50% de suelo desnudo.

El objetivo en este caso es analizar la potencial información que aporta el uso de diferentes índices espectrales y evaluar la información que recuperan sobre las características biofísicas de los ambientes y su variación temporal.

Se utilizará la imagen con reflectancia en superficie de Landsat y vamos a cálcular varios indices. También puede usar otra más que le parezca interesante.

Índice de vegetación normalizado (NDVI)

El índice de vegetación normalizado (NDVI) ayuda a diferenciar la vegetación de otros tipos de cubierta terrestre (artificial) y a determinar su estado general. También permite definir y visualizar las áreas con vegetación en el mapa, así como detectar cambios anormales en el proceso de crecimiento.

vi <- function(img, k, i) {
  bk <- img[[k]]
  bi <- img[[i]]
  vi <- (bk - bi) / (bk + bi)
  return(vi)
}


ndvi <- vi(imagen_stack, 5, 4)
plot(ndvi, col=rev(terrain.colors(10)), main = "Chaco-NDVI in 2022")

Índice de vegetación mejorado (EVI)

El Enhanced vegetation index (EVI) o Índice de Vegetación Mejorado intenta expresar los efectos atmosféricos calculando la diferencia de radiancia entre las bandas del Azul y Rojo y nos permite monitorizar el estado de la vegetación en caso de altas densidades de biomasa.

vi2 <- function(img, m, k, i) {
  bm <- img[[m]]
  bk <- img[[k]]
  bi <- img[[i]]
  vi2 <- 2.5 * ((bm - bk) / (bm + 6* bk - 7.5* bi + 1))
  return(vi2)
}

# For Landsat 8: NIR = 5, red = 4, blue=2.
evi <- vi2(imagen_stack, 5, 4, 2)
  plot(evi, col=rev(terrain.colors(6)), main = "Chaco-EVI in 2022")

par(mfrow=c(1,2))
hist(ndvi, main="NDVI")
hist(evi, main="EVI")

# correlation plot between indices

indices_ndvi<- as.data.frame(ndvi, xy=TRUE, na.rm=TRUE) # convertir a data frame
indices_evi<- as.data.frame(evi, xy=TRUE, na.rm=TRUE) # convertir a data frame
indices <- merge(indices_ndvi,indices_evi,by=c("x","y"))
names(indices) <- c('x','y',"ndvi", "evi")
ggplot(indices, aes(x=evi, y=ndvi)) + 
  geom_point()+
  geom_smooth(method=lm)
## `geom_smooth()` using formula 'y ~ x'

#plot(indices_ndvi$layer,indices_evi$layer)

Índice Diferencial de Agua Normalizado (NDWI)

El NDWI se utiliza como una medida de la cantidad de agua que posee la vegetación o el nivel de saturación de humedad que posee el suelo. Para obtener este índice la combinación de bandas es la siguiente: Landsat 9 (3-6)/(3+6).

wi <- function(img, j, k) {
  bj <- img[[j]]
  bk <- img[[k]]
  wi <-  (bj - bk) / (bj +bk)
  return(wi)
}
ndwi <- wi(imagen_stack, 3,5)
plot(ndwi, col=rev(terrain.colors(10)), main = "Chaco-NDWI in 2022")